Регрессия - метод предсказания (объяснения) одной переменной через другие.
\[y_i = \hat\beta_0 + \hat\beta_1 \times x_i + \epsilon_i,\] * независимые переменные (предикторы, effect) * зависимая переменная (outcome)
В отличие от корреляции, проверяющей связь между двумя переменными, но не направление зависимости, регрессия явно эксплицирует это направление.
library(tidyverse)
library(lme4)
library(lmerTest) # for the model summary
library(skimr) # just for summaries
set.seed(42)
data <- tibble(x = rnorm(150)+1) %>%
mutate(y = 5*x+10+rnorm(100, sd = 2))
mean_y <- mean(data$y)
data %>%
ggplot(aes(x, y)) +
geom_point() +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 0, linetype = 2) +
geom_hline(yintercept = 15.314, color = "blue") +
annotate(geom = "point", x = 0, y = 15.314, size = 4, color = "red") +
scale_x_continuous(breaks = -2:4) +
scale_y_continuous(breaks = c(0:3*10, 15)) +
theme_minimal()
tibble(x = rnorm(150)+1) %>%
mutate(y = 5*x+10+rnorm(100, sd = 2)) %>%
ggplot(aes(x, y)) +
geom_point() +
geom_hline(yintercept = 0, linetype = 2) +
geom_vline(xintercept = 0, linetype = 2) +
geom_smooth(method = "lm", se = FALSE) +
annotate(geom = "point", x = 0, y = 10, size = 4, color = "red") +
annotate(geom = "label", x = -0.5, y = 12, size = 5, color = "red", label = "intercept") +
annotate(geom = "label", x = 2, y = 12.5, size = 5, color = "red", label = "slope") +
scale_x_continuous(breaks = -2:4) +
scale_y_continuous(breaks = c(0:3*10, 15)) +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
pBrackets::grid.brackets(815, 420, 815, 545, lwd=2, col="red")
Когда мы пытаемся научиться предсказывать данные одной переменной \(Y\) при помощи другой переменной \(X\), мы получаем формулу:
\[y_i = \hat\beta_0 + \hat\beta_1 \times x_i + \epsilon_i,\] где
Причем, иногда мы можем один или другой параметр считать равным нулю.
Определите по графику формулу синей прямой.
Задача регрессии — оценить параметры \(\hat\beta_0\) и \(\hat\beta_1\), если нам известны все значения \(x_i\) и \(y_i\) и мы пытаемся минимизировать значния \(\epsilon_i\). В данном конкретном случае, задачу можно решить аналитически и получить следующие формулы:
\[\hat\beta_1 = \frac{(\sum_{i=1}^n x_i\times y_i)-n\times\bar x \times \bar y}{\sum_{i = 1}^n(x_i-\bar x)^2}\]
\[\hat\beta_0 = \bar y - \hat\beta_1\times\bar x\]
При этом, вне зависимости от статистической школы, у регрессии есть свои ограничения на применение:
количество слов и в рассказах М. Зощенко в зависимости от длины рассказа:
zo <- read_tsv("https://github.com/agricolamz/DS_for_DH/raw/master/data/tidy_zoshenko.csv")
zo %>%
filter(word == "и") %>%
distinct() %>%
ggplot(aes(n_words, n))+
geom_point()+
labs(x = "количество слов в рассказе",
y = "количество и")
Избавимся от выбросов и добавим регрессионную линию при помощи
функции geom_smooth():
zo %>%
filter(word == "и",
n_words < 1500) %>%
distinct() ->
zo_filtered
zo_filtered %>%
ggplot(aes(n_words, n))+
geom_point()+
geom_smooth(method = "lm", se = FALSE)+
labs(x = "количество слов в рассказе",
y = "количество и")
Чтобы получить формулу этой линии, нужно запустить функцию, которая оценивает линейную регрессию:
fit0 <- lm(n ~ n_words, data = zo)
fit0
##
## Call:
## lm(formula = n ~ n_words, data = zo)
##
## Coefficients:
## (Intercept) n_words
## 1.259006 0.006204
fit <- lm(n~n_words, data = zo_filtered)
fit
##
## Call:
## lm(formula = n ~ n_words, data = zo_filtered)
##
## Coefficients:
## (Intercept) n_words
## -1.47184 0.04405
\[n_{fit0} = 1.259006 + 0.006204 \times n\_words\]
\[n = -1.47184 + 0.04405 \times n\_words\]
Более подробная информация - в функции summary():
summary(fit0)
##
## Call:
## lm(formula = n ~ n_words, data = zo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.241 -4.974 -2.995 0.385 133.759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.259e+00 1.071e-01 11.76 <2e-16 ***
## n_words 6.204e-03 8.304e-05 74.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.27 on 46656 degrees of freedom
## Multiple R-squared: 0.1068, Adjusted R-squared: 0.1068
## F-statistic: 5581 on 1 and 46656 DF, p-value: < 2.2e-16
summary(fit)
##
## Call:
## lm(formula = n ~ n_words, data = zo_filtered)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.6830 -4.3835 0.8986 4.6486 19.6413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.471840 2.467149 -0.597 0.553
## n_words 0.044049 0.003666 12.015 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.945 on 64 degrees of freedom
## Multiple R-squared: 0.6928, Adjusted R-squared: 0.688
## F-statistic: 144.4 on 1 and 64 DF, p-value: < 2.2e-16
В разделе Coefficients содержится информация о наших
коэффициентах:
Estimate – полученная оценка коэффициентов;Std. Error – стандартная ошибка среднего;t value – \(t\)-статистика, полученная при проведении
одновыборочного \(t\)-теста,
сравнивающего данный коэфициент с 0;Pr(>|t|) – полученное \(p\)-значение;Multiple R-squared и Adjusted R-squared —
одна из оценок модели, показывает связь между переменными. Без поправок
совпадает с квадратом коэффициента корреляции Пирсона:cor(zo_filtered$n_words, zo_filtered$n)^2
## [1] 0.6928376
F-statistic — \(F\)-статистика полученная при проведении
теста, проверяющего, не являются ли хотя бы один из коэффицинтов
статистически значимо отличается от нуля. Совпадает с результатами
дисперсионного анализа (ANOVA).Сколько будет “и” в рассказе Зощенко длиной 1000 слов? (на самом деле, такого рассказа нет. Это предсказание)
predict(fit, tibble(n_words = 1000))
## 1
## 42.57715
Линейная регрессия предполагает нормальность распределения остатков. Когда связь не линейна, то остатки тоже будут распределены не нормально.
Можно смотреть на первый график используя функцию plot()
— график остатков.
plot(fit)
Интерпретаций этого графика достаточно много (см. статью
про это). * Residuals vs Fitted: тренд в стандартизованных остатках
(красная линия) должен идти примерно вдоль пунктирной линии y=0
(нормальность распределения остатков).
* QQplot
* Scale-location plot: показывает гомоскидастичность остатков
(постоянство распределения) - красная линия должна идти примерно
горизонтально.
* Residuals vs. Leverage: показывает точки наблюдения, наиболее сильно
влияющие на коэффициенты модели. Если какая-либо точка вываливается за
границы красных пунктирных линий (Cook’s distance), то ее влияние
чрезмерно.
Можно смотреть на QQplot:
tibble(res = m1$residuals) %>%
ggplot(aes(res))+
geom_histogram(aes(y = ..density..))+
stat_function(fun = dnorm, args = list(mean = 0, sd = sd(m1$residuals)), color = "red")
qqnorm(m1$residuals)
qqline(m1$residuals)
tibble(res = m2$residuals) %>%
ggplot(aes(res))+
geom_histogram(aes(y = ..density..))+
stat_function(fun = dnorm, args = list(mean = 0, sd = sd(m2$residuals)), color = "red")
qqnorm(m2$residuals)
qqline(m2$residuals)
tibble(res = m3$residuals) %>%
ggplot(aes(res))+
geom_histogram(aes(y = ..density..))+
stat_function(fun = dnorm, args = list(mean = 0, sd = sd(m3$residuals)), color = "red")
qqnorm(m3$residuals)
qqline(m3$residuals)
Распределение остатков непостоянно (т.е. не гомоскидастично): остатки становятся больше с ростом значений предиктора (или наоборот).
ldt %>%
ggplot(aes(Mean_RT, Freq))+
geom_point()+
theme_bw()
Тоже решается преобразованием данных.
Линейная связь между некоторыми предикторами в модели.
car::vif()
Наблюдения должны быть независимы. В ином случае нужно использовать модель со смешанными эффектами.
\[y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy_chekhov} + \epsilon_i,\] \[y_i = \hat\beta_0 + \hat\beta_1 \times 1 + \epsilon_i,\]
dummy_chekhov = 1, если это рассказ Чехова dummy_chekhov = 0, если это рассказ Зощенко (дефолтный)
\[ mean(y_z) = \hat\beta_0 \]
\[ mean(y_{ch}) \]
Давайте рассмотрим простой пример с рассказами Чехова и Зощенко, которые мы рассматривали в прошлом разделе. Мы будем анализировать логарифм доли слов деньги в текстах:
chekhov <- read_tsv("https://github.com/agricolamz/DS_for_DH/raw/master/data/tidy_chekhov.tsv")
zoshenko <- read_tsv("https://github.com/agricolamz/DS_for_DH/raw/master/data/tidy_zoshenko.csv")
chekhov$author <- "Чехов"
zoshenko$author <- "Зощенко"
chekhov_zoshenko <-
chekhov %>%
bind_rows(zoshenko) %>%
filter(str_detect(word, "деньг")) %>%
group_by(author, titles, n_words) %>%
summarise(n = sum(n)) %>%
mutate(log_ratio = log(n/n_words)) %>%
ungroup()
Визуализация выглядит так:
chekhov_zoshenko %>%
group_by(author) %>%
mutate(mean = mean(log_ratio)) %>%
ggplot(aes(author, log_ratio))+
geom_violin()+
geom_hline(aes(yintercept = mean), linetype = 3)+
geom_point(aes(y = mean), color = "red", size = 5)+
scale_y_continuous(breaks = c(-7, -5, -3, -6.34))
Красной точкой обозначены средние значения, так что мы видим, что между двумя писателями есть разница, но является ли она статистически значимой? В прошлом разделе мы рассмотрели, что в таком случае можно сделать t-test:
t.test(log_ratio~author,
data = chekhov_zoshenko,
var.equal =TRUE) # поправка Уэлча отключена
##
## Two Sample t-test
##
## data: log_ratio by author
## t = 5.6871, df = 125, p-value = 8.665e-08
## alternative hypothesis: true difference in means between group Зощенко and group Чехов is not equal to 0
## 95 percent confidence interval:
## 0.8606107 1.7793181
## sample estimates:
## mean in group Зощенко mean in group Чехов
## -5.021262 -6.341226
Разница между группами является статистически значимой (t(125) = 5.6871, p-value = 8.665e-08).
Для того, чтобы запустить регрессию на категориальных данных категориальная переменная автоматически разбивается на группу бинарных dummy-переменных:
tibble(author = c("Чехов", "Зощенко"),
dummy_chekhov = c(1, 0),
dummy_zoshenko = c(0, 1))
## # A tibble: 2 × 3
## author dummy_chekhov dummy_zoshenko
## <chr> <dbl> <dbl>
## 1 Чехов 1 0
## 2 Зощенко 0 1
Дальше для регрессионного анализа выкидывают одну из переменных, так как иначе модель не сойдется (dummy-переменных всегда n-1, где n — количество категорий в переменной).
tibble(author = c("Чехов", "Зощенко"),
dummy_chekhov = c(1, 0))
## # A tibble: 2 × 2
## author dummy_chekhov
## <chr> <dbl>
## 1 Чехов 1
## 2 Зощенко 0
Если переменная dummy_chekhov принимает значение 1,
значит речь о рассказе Чехова, а если принимает значение 0, то о
рассказе Зощенко. Если вставить нашу переменную в регрессионную формулу
получится следующее:
\[y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy_chekhov} + \epsilon_i,\]
Так как dummy_chekhov принимает либо значение 1, либо
значение 0, то получается, что модель предсказывает лишь два
значения:
\[y_i = \left\{\begin{array}{ll}\hat\beta_0 + \hat\beta_1 \times 1 + \epsilon_i = \hat\beta_0 + \hat\beta_1 + \epsilon_i\text{, если рассказ Чехова}\\ \hat\beta_0 + \hat\beta_1 \times 0 + \epsilon_i = \hat\beta_0 + \epsilon_i\text{, если рассказ Зощенко} \end{array}\right.\]
Таким образом, получается, что свободный член \(\beta_0\) и угловой коэффициент \(\beta_1\) в регрессии с категориальной переменной получает другую интерпретацию. Одно из значений переменной кодируется при помощи \(\beta_0\), а сумма коэффициентов \(\beta_0+\beta_1\) дают другое значение переменной. Так что \(\beta_1\) – это разница между оценками двух значений переменной.
Запустим регрессию на этих же данных:
fit2 <- lm(log_ratio~author, data = chekhov_zoshenko)
summary(fit2)
##
## Call:
## lm(formula = log_ratio ~ author, data = chekhov_zoshenko)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8652 -0.6105 -0.0607 0.6546 3.2398
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.0213 0.2120 -23.680 < 2e-16 ***
## authorЧехов -1.3200 0.2321 -5.687 8.67e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9717 on 125 degrees of freedom
## Multiple R-squared: 0.2056, Adjusted R-squared: 0.1992
## F-statistic: 32.34 on 1 and 125 DF, p-value: 8.665e-08
Во-первых, стоит обратить внимание на то, что R сам преобразовал нашу
категориальную переменную в dummy-переменную authorЧехов.
Во-вторых, можно заметить, что значения t-статистики и p-value совпадают
с результатами полученными нами в t-тесте выше. Статистически значимый
коэффициент при аргументе authorЧехов следует
интерпретировать как разницу средних между логарифмом долей в рассказах
Чехова и Зощенко.
Дамми-кодирование с более чем двумя уровнями: \[
y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy_chekhov} +
\hat\beta_2 \times \text{dummy_tolstoy} + \epsilon_i,
\] Зощенко Чехова Толстой
dummy_zoschenko 1 0 0
dummy_chekhov 0 1 0
dummy_tolstoy 0 0 1
1 категориальный предиктор (author):
Multiple R-squared: 0.2056, Adjusted R-squared: 0.1992
1 категориальный (author) + 1 непрерывный предиктор (n_word):
Multiple R-squared: 0.4276, Adjusted R-squared: 0.4184
fit3 <- lm(log_ratio ~ author + n_words, data = chekhov_zoshenko)
summary(fit3)
##
## Call:
## lm(formula = log_ratio ~ author + n_words, data = chekhov_zoshenko)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5996 -0.6427 -0.0116 0.5475 3.2360
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.487e+00 1.964e-01 -22.841 < 2e-16 ***
## authorЧехов -1.088e+00 2.006e-01 -5.422 2.95e-07 ***
## n_words -6.873e-04 9.909e-05 -6.936 1.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8281 on 124 degrees of freedom
## Multiple R-squared: 0.4276, Adjusted R-squared: 0.4184
## F-statistic: 46.32 on 2 and 124 DF, p-value: 9.465e-16
Метрика мультиколлинеарности предикторов
car::vif(fit3)
## author n_words
## 1.028633 1.028633
#VIF: 1 / 1-R2
\[R^2 = 1 − Unexplained\_Variation /
Total\_Variation\] Примечание: В случае, если один из предикторов
- непрерывный, а другой - категориальный (фактор), проверяем количество
уровней в факторе. Если уровней всего два, как в данном случае, то
vif() интерпретируется обычным образом. Если уровне более
двух, то выдаются метрики GVIF and aGSIF (generalized VIF). aGSIF >
1.6 говорит о средней инфляции, aGSIF > 2.2 или aGSIF > 3.2
свидетельствует о большой инфляции.
# which components are in lm?
names(fit3)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "contrasts" "xlevels" "call" "terms"
## [13] "model"
Коэффициенты модели:
fit3$coefficients
## (Intercept) authorЧехов n_words
## -4.4869205145 -1.0878105849 -0.0006872767
summary(fit3)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.4869205145 1.964436e-01 -22.840753 2.951741e-46
## authorЧехов -1.0878105849 2.006117e-01 -5.422467 2.951167e-07
## n_words -0.0006872767 9.908636e-05 -6.936138 1.971634e-10
Предсказанные значения зависимой переменной:
head(fit3$fitted.values)
## 1 2 3 4 5 6
## -4.893788 -4.703413 -5.000316 -4.990007 -6.315764 -4.839493
Ошибки модели (остатки, unexplained residuals):
head(fit3$residuals)
## 1 2 3 4 5 6
## -1.48971833 0.33713439 0.94920037 0.59145111 -1.57069351 -0.01448803
Сумма предсказанных значений и ошибок дает значение зависимой переменной:
chekhov_zoshenko %>%
select(log_ratio) %>%
mutate(fit3.predicted = fit3$fitted.values,
fit3.residuals = fit3$residuals,
`predicted+residuals` = fit3$fitted.values + fit3$residuals) %>%
head()
## # A tibble: 6 × 4
## log_ratio fit3.predicted fit3.residuals `predicted+residuals`
## <dbl> <dbl> <dbl> <dbl>
## 1 -6.38 -4.89 -1.49 -6.38
## 2 -4.37 -4.70 0.337 -4.37
## 3 -4.05 -5.00 0.949 -4.05
## 4 -4.40 -4.99 0.591 -4.40
## 5 -7.89 -6.32 -1.57 -7.89
## 6 -4.85 -4.84 -0.0145 -4.85
Дисперсия значений зависимой переменной делится на объясненную моделью и необъясненную. Полная дисперсия (total sum of squares, TSS) может быть подсчитана как сумма квадратов разниц со средним.
tss <- sum((chekhov_zoshenko$log_ratio-mean(chekhov_zoshenko$log_ratio))^2)
Необъясненная дисперсия - сумма квадратов ошибок:
rss <- sum(fit3$residuals^2)
Посчитаем \(R^2\) – долю объясненной дисперсии:
1 - rss/tss
## [1] 0.4276272
#summary(fit3)$r.squared
Модель объясняет 43% дисперсии - много это или мало?
Сравним R2 в моделях fit3 и fit2:
summary(fit2)$r.squared
## [1] 0.2055557
Функция aov() может использоваться для представления выдачи модели в формате выдачи дисперсионного анализа.
aov(fit3)
## Call:
## aov(formula = fit3)
##
## Terms:
## author n_words Residuals
## Sum of Squares 30.53837 32.99204 85.03455
## Deg. of Freedom 1 1 124
##
## Residual standard error: 0.8281078
## Estimated effects may be unbalanced
summary(aov(fit3))
## Df Sum Sq Mean Sq F value Pr(>F)
## author 1 30.54 30.54 44.53 7.45e-10 ***
## n_words 1 32.99 32.99 48.11 1.97e-10 ***
## Residuals 124 85.03 0.69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Хорошая модель должна обладать свойством парсимонией – не только хорошо обобщать данные, но быть настолько сложной, насколько это необходимо (т. е. не иметь больше предикторов, чем это минимально необходимо). Более сложная модель должна доказать, что она значимо лучше объясняет дисперсию. Если это не так, стоит использовать более простую модель.
anova(fit2,fit3)
## Analysis of Variance Table
##
## Model 1: log_ratio ~ author
## Model 2: log_ratio ~ author + n_words
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 125 118.027
## 2 124 85.035 1 32.992 48.11 1.972e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p-value показывает статистическую значимость добавления новых предикторов. Степень свободы (Df) 1 говорит о том, что две модели различаются на одну переменную.
\[y_i = \hat\beta_0 + \hat\beta_1 \times x_{1i}+ \dots+ \hat\beta_n \times x_{ni} + \epsilon_i,\]
В такой регресии предикторы могут быть как числовыми, так и категориальными (со всеми вытекающими последствиями, которые мы обсудили в предудщем разделе). Такую регрессию чаще всего сложно визуализировать, так как в одну регрессионную линию вкладываются сразу несколько переменных.
Попробуем предсказать длину лепестка на основе длины чашелистик и вида ириса:
iris %>%
ggplot(aes(Sepal.Length, Petal.Length, color = Species))+
geom_point()
Запустим регрессию:
fit3 <- lm(Petal.Length ~ Sepal.Length+ Species, data = iris)
summary(fit3)
##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76390 -0.17875 0.00716 0.17461 0.79954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.70234 0.23013 -7.397 1.01e-11 ***
## Sepal.Length 0.63211 0.04527 13.962 < 2e-16 ***
## Speciesversicolor 2.21014 0.07047 31.362 < 2e-16 ***
## Speciesvirginica 3.09000 0.09123 33.870 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2826 on 146 degrees of freedom
## Multiple R-squared: 0.9749, Adjusted R-squared: 0.9744
## F-statistic: 1890 on 3 and 146 DF, p-value: < 2.2e-16
Все предикторы статистически значимы. Давайте посмотрим предсказания модели для всех наблюдений:
iris %>%
mutate(prediction = predict(fit3)) %>%
ggplot(aes(Sepal.Length, prediction, color = Species))+
geom_point()
Всегда имеет смысл визуализировать, что нам говорит наша модель. Если
использовать пакет ggeffects (или предшествовавший ему
пакет effects), это можно сделать не сильно задумываясь,
как это делать:
library(ggeffects)
plot(ggpredict(fit3, terms = c("Sepal.Length", "Species")))
Как видно из графиков, наша модель имеет одинаковые угловые коэффициенты (slope) для каждого из видов ириса и разные свободные члены (intercept).
summary(fit3)
##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76390 -0.17875 0.00716 0.17461 0.79954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.70234 0.23013 -7.397 1.01e-11 ***
## Sepal.Length 0.63211 0.04527 13.962 < 2e-16 ***
## Speciesversicolor 2.21014 0.07047 31.362 < 2e-16 ***
## Speciesvirginica 3.09000 0.09123 33.870 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2826 on 146 degrees of freedom
## Multiple R-squared: 0.9749, Adjusted R-squared: 0.9744
## F-statistic: 1890 on 3 and 146 DF, p-value: < 2.2e-16
\[y_i = \left\{\begin{array}{ll} -1.70234 + 0.63211 \times \text{Sepal.Length} + \epsilon_i\text{, если вид setosa}\\ -1.70234 + 2.2101 + 0.63211 \times \text{Sepal.Length} + \epsilon_i\text{, если вид versicolor} \\ -1.70234 + 3.09 + 0.63211 \times \text{Sepal.Length} + \epsilon_i\text{, если вид virginica} \end{array}\right.\]
Давайте восползуемся данными из пакета Rling
Натальи Левшиной. В датасете 100 произвольно выбранных слов из
проекта English Lexicon Project (Balota et al. 2007), их длина, среднее
время реакции и частота в корпусе.
ldt <- read_csv("https://goo.gl/ToxfU6")
ldt
## # A tibble: 100 × 3
## Length Freq Mean_RT
## <dbl> <dbl> <dbl>
## 1 8 131 819.
## 2 10 82 978.
## 3 7 0 908.
## 4 6 592 766.
## 5 12 2 1125.
## 6 12 9 948.
## 7 3 14013 642.
## 8 11 15 817.
## 9 11 48 931.
## 10 5 290 654.
## # … with 90 more rows
Давайте посмотрим на простой график:
ldt %>%
ggplot(aes(Mean_RT, Freq))+
geom_point()+
theme_bw()
Регрессия на таких данных будет супер неиформативна:
ldt %>%
ggplot(aes(Mean_RT, Freq))+
geom_point()+
geom_smooth(method = "lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
m1 <- summary(lm(Mean_RT~Freq, data = ldt))
m1
##
## Call:
## lm(formula = Mean_RT ~ Freq, data = ldt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -224.93 -85.42 -30.52 81.90 632.66
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 826.998242 15.229783 54.301 < 2e-16 ***
## Freq -0.005595 0.001486 -3.765 0.000284 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 143.9 on 98 degrees of freedom
## Multiple R-squared: 0.1264, Adjusted R-squared: 0.1174
## F-statistic: 14.17 on 1 and 98 DF, p-value: 0.0002843
ldt %>%
ggplot(aes(Mean_RT, log(Freq)))+
geom_point()+
geom_smooth(method = "lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).
ldt %>%
ggplot(aes(Mean_RT, log(Freq+1)))+
geom_point()+
geom_smooth(method = "lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
m2 <- summary(lm(Mean_RT~log(Freq+1), data = ldt))
m2
##
## Call:
## lm(formula = Mean_RT ~ log(Freq + 1), data = ldt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -242.36 -76.66 -17.49 48.64 630.49
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1001.60 29.79 33.627 < 2e-16 ***
## log(Freq + 1) -34.03 4.76 -7.149 1.58e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 124.8 on 98 degrees of freedom
## Multiple R-squared: 0.3428, Adjusted R-squared: 0.3361
## F-statistic: 51.11 on 1 and 98 DF, p-value: 1.576e-10
m1$adj.r.squared
## [1] 0.1174405
m2$adj.r.squared
## [1] 0.336078
Отлогорифмировать можно и другую переменную.
ldt %>%
ggplot(aes(log(Mean_RT), log(Freq + 1)))+
geom_point()+
geom_smooth(method = "lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
m3 <- summary(lm(log(Mean_RT)~log(Freq+1), data = ldt))
m1$adj.r.squared
## [1] 0.1174405
m2$adj.r.squared
## [1] 0.336078
m3$adj.r.squared
## [1] 0.3838649
Как интерпретировать полученную регрессию с двумя отлогорифмированными значениями?
В обычной линейной регресии мы узнаем отношения между \(x\) и \(y\): \[y_i = \beta_0+\beta_1\times x_i\]
Как изменится \(y_j\), если мы увеличем \(x_i + 1 = x_j\)? \[y_j = \beta_0+\beta_1\times x_j\]
\[y_j - y_i = \beta_0+\beta_1\times x_j - (\beta_0+\beta_1\times x_i) = \beta_1(x_j - x_i)\]
Т. е. \(y\) увеличится на \(\beta_1\) , если \(x\) увеличится на 1. Что же будет с логарифмированными переменными? Как изменится \(y_j\), если мы увеличем \(x_i + 1 = x_j\)?
\[\log(y_j) - \log(y_i) = \beta_1\times (\log(x_j) - \log(x_i))\]
\[\log\left(\frac{y_j}{y_i}\right) = \beta_1\times \log\left(\frac{x_j}{x_i}\right) = \log\left(\left(\frac{x_j}{x_i}\right) ^ {\beta_1}\right)\]
\[\frac{y_j}{y_i}= \left(\frac{x_j}{x_i}\right) ^ {\beta_1}\]
Т. е. \(y\) увеличится на \(\beta_1\) процентов, если \(x\) увеличится на 1 процент.
Логарифмирование — не единственный вид траснформации:
shiny::runGitHub("agricolamz/tukey_transform")
В датасет собрана частотность разных лемм на основании корпуса НКРЯ [Ляшевская, Шаров 2009]. Известно, что частотность слова связана с рангом слова (см. закон Ципфа). Постройте переменную ранга и визуализируйте связь ранга и логорифма частотности с разбивкой по частям речи. Какие части речи так и не приобрели после трансформации “приемлемую” линейную форму?
\[\sigma^2_X = \frac{\sum_{i = 1}^n(x_i - \bar{x})^2}{n - 1},\]
где
Представим, что у нас есть следующие данные:
set.seed(42)
df <- tibble(x = sort(rnorm(20, mean = 50, sd = 10)),
y = seq_along(x))
df %>%
ggplot(aes(x)) +
geom_point(y = 0) +
ggrepel::geom_text_repel(aes(label = y), y = 0) +
labs(x = "значение наблюдений x")
Дисперсия - сумма квадратов расстояний от каждой точки до среднего выборки (пунктирная линия) разделенное на n-1.
Распределения могут иметь одинаковое среднее, но разную дисперсию:
x <- rnorm(20, mean = 50, sd = 10)
var(x)
## [1] 107.8731
sd(x) # sqrt(var(x))
## [1] 10.3862
\[cov(X, Y) = \frac{\sum_{i = 1}^n(x_i - \bar{x})(y_i-\bar{y})}{n - 1},\] \[\rho_{X,Y} = \frac{cov(X, Y)}{\sigma_X\times\sigma_Y} = \frac{1}{n-1}\times\sum_{i = 1}^n\left(\frac{x_i-\bar{x}}{\sigma_X}\times\frac{y_i-\bar{y}}{\sigma_Y}\right),\]
для коэффициента корреляции Пирсона, где
Для таких данных:
Ковариацию можно представить графически следующим образом:
Положительная ковариация - много красных квадратов, отрицательня ковариация - много синих.
Коэффициент корреляции Пирсона можно представить как среднее произведение \(z\)-нормализованных значений двух переменных.
Что можно сказать про корреляцию в следующих данных:
В открытом доступе можно найти игры “Угадай корреляцию” здесь или здесь.